Ground-state properties of hard core bosons in one-dimensional harmonic traps 
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density matrix is defined as 
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The one-particle density matrices for liard core bosons in a one-dimensional harmonic trap are 
computed numerically for systems with up to 160 bosons. Diagonalization of the density matrix 
shows that the many-body ground state is not Bose-Einstein condensed. The ground state occu- 
pation, the amplitude of the lowest natural orbital, and the zero momentum peak height scale as 
powers of the particle number, and the corresponding exponents are related to each other. Close to 
its diagonal, the density matrix for hard core bosons is similar to the one of noninteracting fermions. 
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Bose-Einstein condensation (BEC) in dilute atomic va- 
pors is a very attractive field of contemporary scientific 
research. Of particular interest are one-dimensional sys- 
tems 1^ which are experimentally reaUzed in very elon- 
gated three-dimensional traps The one-dimensional 
Bose gas displays a rich behavior ranging from the 
weakly interacting regime ||^ to the experimentally re- 
alized Thomas-Fermi regime to a gas of impenetrable 
hard core bosons 0, ||, ^ at sufficiently low densities and 
large s-wave scattering length. This latter case of hard 
core interactions is particularly interesting since the ex- 
act many-body ground-state wave function [See Eq. (Q)] 
is known due to a boson-fcrmion mapping M]. In spite 
of this knowledge, it is nontrivial to extract the ground- 
state properties from this wave function. The main diffi- 
culty to overcome is the computation of the one-particle 
density matrix for sufhciently large systems. Girardeau 
et al. Q used Monte Carlo integration techniques to 
compute the density matrix for small systems with up 
to iV = 10 bosons. More accurate methods employed by 
Lapeyre et al. Q are restricted to up to = 8 particles. 
Analytical results were derived by Kolomeisky et al. |^ 
for the diagonal of the density matrix and by Minguzzi 
et al. for the tails of the momentum distribution. It 
is the purpose of this work to compute the one particle 
density matrix for large systems containing more than a 
hundred particles and to present analytical results and 
scaling relations as well. 

The ground-state wave function for N hard core bosons 
is given by 

N 



PB{x,y) ^ ^ J dzi...dzN-iip*Bizi,...,ZN-i,x) 

xipBizi,...,ZN-i,y) (3) 

and requires the integration over — 1 variables. Let us 
rewrite the one-particle density matrix as 



PB{.x,y) 



,i/2(^_i),°-p|-2(- +r) 

X / d^N-i W \zk - x\\zk - y\. (4) 
fe=i 



Here we introduce the integration measure 



k=l 



l<i<j<N-l 



xl;B{.z,,...,ZM)^Cl\{e-^-'^ [] \z,-z,\ (1) 



k=l l<i<j<N 

where the normalization constant is 



(2) 



Note that the wave function (|^) is simply the absolute 
value of the ground-state wave function ijjp for N nonin- 
teracting fermions |Q, i.e. ipB = I'i/'fI- The one-particle 



which is identical to the joint probability density for 
eigenvalues zi, . . . , z^r-i of {N — l)-dimensional random 
matrices drawn from the Gaussian Unitary Ensemble 
(GUE) |ll|,|l|. 

This observation greatly facilitates the Monte Carlo in- 
tegration since it is easy to generate independent configu- 
rations zi, . . . , zn-i that are distributed according to the 
measure (^. In practice, one draws a (A^— l)-dimensional 
random matrix from the GUE, i.e. the matrix is complex 
Hermitian, and its matrix elements are Gaussian random 
variables with zero mean and variance 1/2. Upon diago- 
nalization, one obtains eigenvalues zi, . . . , zn-i that are 
distributed according to the measure (||). The integrand 
rife^Ji^ |zfe — x\\zk — y\ is then evaluated for this configu- 
ration, and the procedure is repeated many times. This 
procedure requires an effort of 0{N^) for each configura- 
tion. The computation of the density matrix for N — 10 
bosons takes about one hour on a PC and requires 10^ 
configurations. While this is already a considerable im- 
provement over previous integration techniques ^, ||] , it 
is not sufficient when much larger systems are considered. 
Following Mehta jl^, we express the measure (||) in 



2 



terms of harmonic oscillator wave functions ipn {z) as 
dzi . . . dzpf. ' ^2 



5.59 



(iV-1) 



^yr^ (det [<Pfc-i(2;i)]fc,i=i,...,iv-i) 



(6) 

The integration can be performed since the [N — 1)- 
dimensional integral factorizes 



dzx ■ ■ .dzN~i (det [(pk^i{zi)\zi - x\]|^ l^-^^ J^_-^ 

X (dct[(p,^iizj)\Zj -2;|]ij = i,._Ar_i) 

= (Ar~l)!det [Bm,nix,y)] 



m,n=0,...,N-2 ' 



(7) 



Here we introduce the [N — l)-dimensional square matrix 
with elements 



Bm.,nix,y) = / dz\z - X\\z - y\iprniz)ipn{z). (8) 



We may thus express the one-particle density matrix of 
N bosons as the determinant of a {N — l)-dimensional 
matrix 



Pb{x,v) 



7ri/2(iV_ 1)1 



det [Bra,n{x,v)\ 



m.n=0,...,N-2 ' 



(9) 

This form of the density matrix and its connection to 
the GUE was previously discussed by Forrester et al. . 
The density matrix pB{x,y) may now be computed nu- 
merically. For a numerical computation of Hermite poly- 
nomials, see, e.g., Refs.jlJ, The total effort scales 
like 0{N^), and the calculation for a system of iV = 100 
bosons requires about one day on a PC. We computed 
the density matrix for systems containing up to iV = 160 
particles. Figures |l| and || show density plots of the den- 
sity matrix pB{x,y) for iV = 10 and N = 100 particles, 
respectively. The salient features are strong intensities 
close to the diagonal x — y and small off-diagonal contri- 
butions which rapidly drop to zero at \x\, \y \ « R, where 



R = {2Nf/'^ 



(10) 



is the radius of the density psix, x) [gj or Wigner's semi- 
circle [|ll|. Note that the width of the peak along the 
diagonal decreases with increasing particle number. 

Diagonalization of the density matrix yields the nat- 
ural orbitals 4>j{x) and occupation numbers > ni > 
712 > . . . > fulfilling / dypB{x,y)(t>j{y) = nj(j)j{x), i.e. 



PB{x,y) = ^nj(t)j{x)<t>j{y). 



(11) 



The many-body ground-state is Bose-Einstein condensed 
provided the density matrix exhibits a macroscopic eigen- 
value no oc N; the corresponding natural orbital (/>o(x) 
is then the wave function of the BEC. Our data show 
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FIG. 1: One-particle density matrix psip^^y) for A'^ = 10 
bosons. Lightest gray indicates almost zero amplitude, while 
black indicates the maximal amplitude « {2NY^'^ /-k. The 
width of the diagonal peak is proportional to N~^^^. 



that no is not simply a power of N for the limited range 
of particle numbers considered in this work. Previous 
studies on small systems suggested that no cx N°' with 
a « 0.59 1^. While this value is a good fit for systems 
with up to TV = 10 particles, we find that the exponent 
a actually decreases with increasing particle number N, 
reaching a ~ 0.53 for the largest particle numbers consid- 
ered in this work. Below we will present simple scaling 
arguments which suggest that a = 0.5 is the expected 
behavior for large particle number N. This behavior of 
harmonically trapped hard core bosons is similar to the 
uniform system of hard core bosons p^ ] . 

Figure ^ shows the scaled natural orbitals R^^'^<po{x) 
for N — 10, 40, 160 bosons. Note that the scaled nat- 
ural orbitals approach an A^-independent function that 
depends only on the scaled variable x/R as N increases. 
This function is nonzero only for ^ 1 and varies 

smoothly over this interval. We find that the maxi- 
mum amplitude of the lowest natural orbital scales like 
(/>o(0) cx with P = -0.25. This is expected since the 
natural orbital 0o is normalized and is supported on a 
domain that scales like R oc N^^^. 

The momentum distribution 



n(fc) = 



1 

2^ 



dxdypB{x,y)e 



(12) 



is of particular interest. Plots of the normalized mo- 
mentum distribution n(k)/N are shown in Figure I for 
N — 10,40,160. The distributions have a pronounced 
peak at zero momentum and long tails. Minguzzi et al. 
pO| showed that the tails of the momentum distribution 
decay as n(fc) cx for large momenta k. We find that 
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FIG. 2: Same as Figure ^ but for N = 100 particles. 



the peak height n(0) is proportional to the particle num- 
ber, n{0) cx iV^ with 7 = 1.0 (see the inset of Figure |^). 
Thus, the system of hard core bosons mimics the macro- 
scopic occupation of a momentum zero state and in this 
aspect resembles a uniform and noninteracting Bose sys- 
tem. 

This finding is particularly interesting, because it al- 
lows us to predict the iV-dependence of the ground-state 
occupation tiq- According to Eq. (|l2|), the momentum 
peak height n(0) is simply the integral over the density 
matrix which we approximate as 



n(0) oc J dxdy pB{x,y) 
diagonal region 

+ 110 j dxdy(j)o{x)(j)oiy). (13) 
off-diagonal domain 

Here we have decomposed the domain of integration into 
the diagonal region over the peaked structure and the off- 
diagonal contribution; the latter may be approximated 
in leading order by the integral over the lowest natural 
orbital (j>o since all other natural orbitals have smaller 
occupation number and their positive and negative am- 
plitudes lead to cancellations when integrated over the 
off-diagonal domain. The contribution from diagonal re- 
gion scales like A^^/^ and can therefore be neglected when 
compared with the momentum peak n(0) oc iV^'°. This 
can be seen as follows. The normalized density along 
the diagonal is proportional to TV, while the width of 
the diagonal peak decreases like N~^^^ with increasing 
particle number (see below). The remaining off-diagonal 
contribution scales like noi?^((/>o(0))^ = noN^-^ since the 
domain of integration is a square whose area scales like 
(X N . Thus, no oc iV" with a = 0.5. Let us summa- 




FIG. 3: Scaled natural orbital R^^^<j}o{x) for N ^ 10 (full 
line), A*' = 40 (dashed-dotted line), and = 160 (dashed 
line) bosons. R = {2NY^'^ sets the length scale in units of 
the oscillator length. The natural orbitals 4>o{^) ^-re nonzero 
on an interval that scales like A'"^''^ and have a maximum 
amplitude proportional to N~''^'^^ . 



rize the analysis of this paragraph into the equation 

7 = l + a + 2/3 (14) 

that relates the scaling of the momentum peak height to 
the ground state occupation and the amplitude of the 
lowest natural orbital. 

The density matrix pB{x,y) can be understood analyt- 
ically close to its diagonal x = y. Kolomeisky et al. j|] 
showed that the diagonal density matrix is identical to 
the density of noninteracting fcrmions in one-dimensional 
harmonic traps, i. e. pb(x,x) = 'Y^n=^ which 
is also identical to the level density of GUE random ma- 
trices ll^. This analogy holds also in leading order as 
one leaves the diagonal. To see this, we write the matrix 
(S) as 



:(a:,y) 



-^m,n(*^i y) — -^rr 



with 



dz{z-x){z-y)ipjn{z)ipn{z), 



(15) 



rn — 1 ) 



/oo 
dz{z - x){z - y)iprn{z)^n{z) 
-00 

+ i (v/m(n +l)<5;\t\ + v/n(m + l)C+\ 

The decomposition of the matrix elements Bm,n into 
the matrix elements Fm,n and a remainder is quite useful. 
The remaining integral in Eq. (15) vanishes on the diag- 



onal X = y and yields corrections proportional to jx — t/|' 
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FIG. 4: Normalized momentum distribution n{k)/N for sys- 
tems oi N = 10 (full line), A'' = 40 (dashed-dotted line), and 

= 160 (dashed line) bosons. The momentum k is given 
in units of the inverse oscillator length. The inset shows the 
zero-momentum peak n(0) as a function of particle number 

in a log-log plot. 



as one leaves the diagonal. The leading corrections pro- 
portional to {x — are already contained in the matrix 
elements F,„_„. Thus, close to the diagonal, the matrix 
B can be approximated by the matrix F. This latter ma- 
trix, however, determines the density matrix pp{x^y) for 
noninteracting fermions in a one-dimensional harmonic 
trap. Indeed, adjusting the above calculation to the case 
of noninteracting fermions yields 

2iV-lg-(i.= +y")/2 
PP^'^^y) = 7rV2(jV^l)! [■^™^»(^'y)L,n=0,...,JV-2 

for the fermionic density matrix. [Note that the more 
familiar expression pp{x,y) = '^n=Q Vn{x)Lpn{y) for the 
same density matrix can be obtained in a direct calcula- 
tion]. Thus, the density matrix for hard core bosons is 
practically identical to the density matrix for noninter- 
acting fermions close to its diagonal. Note, finally, that 
the square of the fermionic density matrix p%{x,y) is 
identical to the two- level cluster function of the GUE . 
This allows us to transfer results obtained for random 
matrices to the case of noninteracting fermions. Close to 
the origin, we thus find for large particle number N ^ 1 
and 5 < 



Pb{5, -5) « pf{5, -6) = -^/2N (l - ^N5^ 

TT V 3 



Thus, the width and height of the prominent peak along 
the diagonal of the bosonic density matrix (see, e.g.. Fig- 
ures ^ |2|) scales like iV~^/^ and like iV^/^, respectively. 
This behavior is confirmed by our numerical computa- 
tions where we take the width at half maximum. 



In summary, we expressed the density matrix of N har- 
monically trapped hard core bosons as a determinant of a 
(A'^— l)-dimensional symmetric matrix and performed nu- 
merical computations on systems containing more than a 
hundred particles. The density matrix is strongly peaked 
along its diagonal but lacks off-diagonal long range or- 
der. Accordingly, the ground-state is not Bose-Einstein 
condensed. The ground-state occupation, the amplitude 
of the lowest natural orbital, and the momentum peak 
height scale as powers of the particle number, and the 
corresponding exponents are related to each other by a 
simple equation. The lowest natural orbital approaches a 
particle number independent function when scaled prop- 
erly. Analytical results show that the density matrices 
of hard core boson systems and noninteracting fermion 
systems are almost identical close to the diagonal. 
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